%%%%%%%%%%%%
%%%%%%%%%%%%
%%%%% GEVOLGDE STAPPEN %%%%%%%
%%%%%%%%%%%%
\section{Gevolgde stappen - Deel 1}
In deze sectie wordt een overzicht gegeven van de gevolgde stappen, met bij elke stap een korte uitleg. 

\subsection{Referentieimplementatie van de drie algoritmes}
Alvorens te beginnen werken aan de FORTRAN implementatie is eerst in MATLAB een referentieimplementatie van de algoritmes \code{taylorexp}, \code{eig} en \code{eigqr} gemaakt, en een testprogramma. Dit was vooral om te kijken of de interpretatie van het algoritme klopte, en om te kijken welke LAPACK routines juist nodig zouden zijn. De LAPACK bibliotheek bleek nodig te  zijn om eigenwaarden en eigenvectoren te berekenen, een QR-ontbinding uit te rekenen, en een stelsel op te lossen. De MATLAB referentieimplementaties zijn later ook gebruikt om tussenresultaten van de ge\"implementeerde algoritmes te vergelijken met de uitkomsten van MATLAB, om fouten te kunnen lokaliseren.
\subsection{Opzoekwerk LAPACK routines, extra argumenten, ...}
De LAPACK routines vermeld in de mail waren xGEQRF en xHSEQR. De methode xHSEQR is uiteindelijk niet gebruikt omdat naast de eigenwaarden ook de eigenvectoren nodig zijn. Als precisie is telkens gekozen voor re\"ele of complexe dubbele precisie (D en Z).
\newline
\newline
Om te werken met command-line argumenten is eerst gekeken naar de methodes \code{iargc()} en \code{getarg()}, maar deze werkten enkel met \code{gfortran}. Uiteindelijk werd gekozen voor \code{command\_argument\_count()} en \code{get\_command\_argument()}.
\newline
\newline
Een volledige lijst met externe routines is te vinden in sectie \ref{sec:extrasoftware}.


\subsection{Structuur van de code schrijven en compileren}
De code voor deel 1 is als volgt gestructureerd\footnote{Een klassediagramma van het volledige project is te vinden in bijlage \ref{sec:klassediagramma}}:
\begin{itemize}
\item Een module \code{methods.f95} die alle algoritmes bevat, plus enkele private subroutines.
\item Een module \code{IO.f95} die het printen en inlezen van matrices regelt.
\item Een wrapper \code{matrixexp.f95} die afhankelijk van de argumenten de verschillende algoritmes aan kan roepen.
\item Een programma \code{test.f95} dat de resultaten voor verschillende testmatrices vergelijkt met de oplossing van MATLAB.
\end{itemize}
Van de methodes in \code{methods.f95} is eerst de interface aangemaakt. Dan zijn de methodes in \code{IO.f95} toegevoegd, en de verwerking van de argumenten in \code{matrixexp.f95}. Daarna werden de testen ge\"implementeerd in \code{test.f95}. Zo kon bij het implementeren van de algoritmes de juiste werking onmiddellijk gecontroleerd worden.
\subsection{Schrijven van een testprogramma}
De gebruikte testmatrices voor deel 1 zijn voor een deel deze aangereikt in de opgave, en voor een deel zelf aangemaakt met MATLAB. Eerst is het gebruik van system calls overwogen om in het testprogramma de wrapper te kunnen gebruiken. Dit was echter platform-afhankelijk, en werkte ook niet even goed met alle compilers. 

De module \code{test} gebruikt dus de module \code{IO} om de testmatrices in te lezen, gebruikt de module \code{methods} om de algoritmes aan te roepen en drukt de resultaten af op de standaarduitvoer.

Op de resultaten voor de verschillende testmatrices wordt verder ingegaan in sectie \ref{sec:resultaten}.
\subsection{Implementatie van de Taylor benadering}
In het algoritme \code{taylorexp} wordt formule \ref{eqn:taylorexp} expliciet toegepast. De implementatie van deze methode gaf geen noemenswaardige problemen, en de testen gaven goede resultaten.
\subsection{Implementatie van het eenvoudige algoritme}
Bij de implementatie van het algoritme \code{eig} doken al wat meer problemen op, vooral omdat de keuze was om alles meteen in complexe dubbele precisie te implementeren. Een van de grootste problemen was het declareren van dubbele precisie complexe getallen. De functie \code{dcmplx} wordt niet door elke compiler ondersteund. Uiteindelijk was de oplossing om de functie \code{cmplx} te gebruiken met een extra kind-argument.

In de testen gaf het algoritme goede resultaten voor symmetrische matrices, maar hele slechte voor willekeurige matrices.De verwachting was toen dat de onnauwkeurigheid van de methode \code{eig} tov de methode \code{eigqr} aan de basis lag van deze problemen.
\subsection{Implementatie van het uitgebreide algoritme}
\label{sec:eigqrimplem}
De implementatie van het algoritme \code{eigqr} nam het meeste tijd in beslag, en hier kwam de vergelijking met tussenresultaten van de MATLAB referentieimplementaties zeker van pas. Een moeilijkheid hierbij was dat de eigenvectoren voor complexe eigenwaarden genormaliseerd nog steeds kunnen verschillen in een factor $\pm i$, waardoor de tussenresultaten ook konden verschillen. Een andere moeilijkheid was dat eerst de routine ZGEPQF3 gebruikt werd voor de QR-ontbinding, maar deze gebruikte een kolompivotering, waardoor de matrix $Q$ gepermuteerde kolommen kon hebben.
\newline
\newline Nadat de fouten uit het algoritme waren weggewerkt gaven de testen juiste resultaten voor de voorbeelden uit de opgave, maar enorm grote resultaten voor sommige willekeurige matrices. Door vergelijking met de MATLAB resultaten bleek uiteindelijk dat de methode ZGEEV wel correcte eigenwaarden gaf, maar dat de eigenvectoren soms slechts drie juiste cijfers hadden. Hierna werd ZGEEV vervangen door DGEEV, die wel de eigenvectoren tot op dubbele precisie correct teruggaf. Er moest dan wel een extra private routine toegevoegd worden om de complexe eigenwaarden en eigenvectoren samen te stellen. 
\newline
\newline
Door de berekening van de eigenvectoren te verbeteren was er echter nog geen  verbetering merkbaar in de resultaten voor niet-symmetrische willekeurige matrices. Dit werd na vergelijken met MATLAB tussenresultaten teruggebracht tot onnauwkeurigheden in de QR-ontbinding met ZGEQRF, die zich snel voortplanten als ze daarna in een stelsel terechtkomen. 
Ter illustratie van deze fout een vergelijking van de $R$ uit de QR-ontbinding van $V$, zoals berekend in MATLAB, en zoals berekend met ZGEQRF, voor exact dezelfde $V$ (tot op 15 cijfers nauwkeurig). De $A$-matrix was asymmetrisch met dimensie 4x4.
Uit tabel 1 kan dan worden afgeleid dat de fout voor de drie elementen linksboven quasi nul is, en dat de fout rechtsonder heel snel toeneemt, met uitschieters tot $10^9$. Dit verschijnsel treedt enkel op bij complexe $V$, indien de matrix $A$ symmetrisch is en $V$ dus re\"eel, is de QR-ontbinding wel exact.
\begin{table}[htdp]
\begin{center}
\begin{tabular}{cccc}
      exact  & ${-15.3525} $ &${ -0.8648}$ &   ${-0.8648}$ \\
      0    &  exact &  ${ 4.9486}$ &  ${ 4.9486}$ \\
      0     & 0   & ${4.1903} $ &  ${9.1184} $\\
      0     & 0     & 0 &  ${8.2966} $
\end{tabular}
\caption{$Log_{10}$ van het verschil tussen $R$ berekend met MATLAB en $R$ berekend met ZGEQRF.}
\end{center}
\label{tbl:rfout}
\end{table}%
\newline

\subsection{Eenvoudige algoritme revisited}
Door bij de implementatie van \code{eigqr} de eigenwaardenontbinding nader te bekijken, moesten de verbeteringen ook aangebracht worden in \code{eig}. Bij deze methode had de verbeterde nauwkeurigheid van de eigenvectoren wel degelijk resultaat, omdat hierbij geen onnauwkeurige QR-ontbinding moest worden uitgevoerd. De fouten lagen nu in de buurt van de machineprecisie.

\section{Gevolgde stappen - Deel 2}
\subsection{Feedback van Deel 1 verwerken}
Als feedback op Deel 1 werden enkele kleine fouten vermeld, die bij het aanvangen van Deel 2 als eerste werk zijn aangepast.
\subsection{Referentieimplementatie van FFTW}
De grootste moeilijkheid bij de algoritmes voor circulante matrices was het berekenen van de Fouriertransformatie. Om te controleren of deze correct verliep werden referentieimplementaties gemaakt van \code{taylorexpc} en \code{eigc}, de versies voor circulante matrices. Deze referentieimplementaties gaven ook de tussenresultaten van de voorwaartse en achterwaartse Fouriertransformatie terug.
\subsection{Uitbreiden van de \code{IO} module}
\label{sec:deel2IO}
De \code{IO} module moest nu ook circulante matrices kunnen inlezen, met een aparte formatting. Er is gekozen om een methode \code{iscirculant} te implenteren, die bepaalt of een bestand een circulante matrix bevat of niet. Oorspronkelijk gebruikt deze methode een \code{REWIND} statement om de file pointer waar FORTRAN begon te lezen terug aan het begin van de file te zetten. Dit gaf echter problemen wanneer de input van de standaard input stream gelezen moest worden. Daarom is gekozen om de file pointer wel een lijn te laten opschuiven, en is een optioneel argument aan de \code{readmatrix} methodes toegevoegd om een eventuele overbodige eerste regel over te kunnen slaan. Dit argument moet op \code{.true.} gezet worden indien bijvoorbeeld geen call naar \code{iscirculant} is gemaakt.  
\newline
\newline
De interfaces \code{printmatrix} en \code{printresult} zijn ook toegevoegd aan de module \code{IO}. Deze interfaces omvatten methodes om debugging informatie in leesbaar formaat af te drukken, en om volle en circulante matrices in het juiste formaat op een stream af te drukken. Een deel van deze methodes stond in deel 1 nog verspreid over de verschillende modules, maar zijn in deel 2 verzameld in in \code{IO}. Een klassediagramma van het volledige project met daarin ook de publieke methodes van elke module is opgenomen in bijlage \ref{sec:klassediagramma}.
\newline
\newline
Tenslotte werd de \code{IO} module ook aangepast om grote matrices te kunnen inlezen wanneer hun rijen verspreid staan over meerdere regels.
\subsection{Extra routines, FFTW, ...}
Om de fouriertransformatie te berekenen is de FFTW bibliotheek gebruikt. Voor deel 2 waren voor de rest geen extra LAPACK routines nodig. Een volledige lijst met externe routines is te vinden in sectie \ref{sec:extrasoftware}.
\subsection{Uitbreiden van de testen}
In deel 1 bestond de testmodule enkel uit het aanroepen van de verschillende methodes voor enkele testmatrices, en het tonen van de absolute fout in vergelijking met de uitkomst berekend met \code{expm} uit MATLAB. 
Voor deel 2 is getracht om zoveel mogelijke stukken code te bereiken met de testen. Ook is de aanmaakt van de testmatrices zo veel mogelijk geautomatiseerd.
\newline
\newline
Om naast de module \code{methods} ook de modules \code{IO} en \code{fftw} te testen is er een hulpmodule \code{auxtest} aangemaakt. Deze module bevat een subroutine \code{iotest} om de verschillende \code{readmatrix} en \code{printmatrix} te testen. Ook is er een subroutine \code{ffttest} om voor matrices van verschillende groottes de voorwaartse en achterwaartse Fouriertransformaties te vergelijken met de uitkomst van MATLAB. Het programma \code{test} is aangepast om naast de testen uit \code{auxtest} ook de helpboodschap te testen.
\newline
\newline
De aanmaak van testmatrices is geautomatiseerd met een matlabscript \code{creatematrices.m}. Hierdoor kan gemakkelijk de minimale en maximale grootte en de norm van elke matrix aangepast worden. Voor de standaardtesten is gekozen om de matrices te normaliseren door te delen door de 2-norm. De gebruikte MATLAB bestanden zijn te vinden onder \code{/Refimplementatie}.
\newline
\newline
In deel 2 is de weergave van de absolute fout vervangen door de relatieve fout, omdat deze een beter beeld geeft van de nauwkeurigheid. Naast deze absolute fout wordt ook voor elke combinatie testmatrix/algoritme aangegeven of de fout lager is dan een bepaalde drempelwaarde zodat snel kan worden nagekeken of alle testen goed zijn verlopen. Ook wordt voor elke uitvoering van het algoritme de timing weergegeven, zodat eventuele problemen met een te lange uitvoeringstijd gemakkelijk opgespoord kunnen worden. Bij het uitvoeren van de experimenten beschreven in sectie \ref{sec:resultaten} werden de opgemeten fouten en timings ook weggeschreven naar een bestand, zodat deze gemakkelijk konden worden ingevoerd in MATLAB. Dit vereenvoudigde het bepalen van de tijdscomplexiteit tov deel 1, waarin het UNIX-commando \code{time} werd gebruikt in een shell-script.
\newline
\newline
Naast het testen van de modules met het programma \code{test}, zijn achteraf ook nog enkele testen ge\"implementeerd in een shell-script \code{testmatrixexp.sh}. Dit shell-script test de uitvoering van het gecompileerde en gelinkte programma \code{matrixexp}. De verschillende argumenten worden gecontroleerd, en er wordt bijvoorbeeld ook nagegaan of 
\newline 
\newline \hspace{12pt}\code{./matrixexp ---algoritm taylor < A.in | ./matrixexp} 
\newline 
\newline een geldige uitdrukking is. Hiermee is het probleem met de \code{REWIND} functie vermeld in \ref{sec:deel2IO} ontdekt. 
\subsection{Extra commentaar}
Met het groeien van de modules werden de interfaces van de verschillende subroutines en functies ook minder overzichtelijk. Daarom is geopteerd voor een meer formele vorm van documentatie, waarbij de verschillende parameters en het resultaat allemaal expliciet vermeld worden. Op deze manier werden ook alle argumenten nog eens ge\"evalueerd, en werden overbodige argumenten en methodes verwijderd of samengevoegd.
\subsection{EIG en EIGQR revisited}
Op het einde van deel 1 was er nog steeds een probleem met de LAPACK routines. Dit probleem was gedeeltelijk opgelost door LAPACK handmatig te compileren met \code{gfortran} op mijn eigen machine. Het programma gaf dan de juiste resultaten maar kon enkel gecompileerd worden met gfortran. Uiteindelijk zijn handmatig gecompileerde LAPACK bibliotheken beschikbaar gesteld voor \code{gfortran}, \code{g95} en \code{ifort}, waarmee het probleem opgelost was. De compiler \code{nagfor} is niet gebruikt bij de testen bij gebrek aan werkende LAPACK-bibliotheek.
\subsection{Implementatie van de Methodes}
Het basisalgoritme voor de drie methodes zoals beschreven in sectie \ref{sec:theoriedeel2} is zeer eenvoudig, en kan op enkele regels neergeschreven worden:
\lstset{language=Fortran,keywordstyle=\color{blue},label=DescriptiveLabel}
\begin{lstlisting}
function expm(Ac) result(RES)
    real(doubleprecision), dimension(:), intent(in) :: Ac
    real(doubleprecision) :: RES(size(Ac,1))
    complex(doubleprecision) :: c(size(Ac,1))
    c=fft(Ac)
    c=exp(c)
    RES=ifft(c)
end function
\end{lstlisting}
De enige aanpassingen die hier nog aan gemaakt zijn zijn de mogelijkheid om debugging informatie af te drukken, en het berekenen van \code{exp(c)} in \code{taylorexpc} met de taylorexpansie.
\subsection{Aanpassen van matrixexp}
In \code{matrixexp} zijn enkele aanpassingen gebeurd om het inlezen van circulante matrices toe te laten, en om het juiste algoritme te selecteren voor elke input.

